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ABSTRACT 


A new quasilinearization method, referred to 
as the Modified Quasilinearization Method, is pro- 
posed for numerically solving the nonlinear two- 
point boundary value problem with an undetermined 
terminal time. This method is an extension to pre- 
viously proposed quasilinearization methods, in that 
the terminal boundary may be specified by a general 
function of the problem variables and time, rather 
than just by specific values of the variables. More- 
over, for variable terminal- time problems the ter- 
minal time determination is made an integral part of 
the iterative method itself. In addition, a scheme is 
formulated and successfully implemented which sig- 
nificantly reduces the sensitivity of the convergence 
characteristics of the method to the required initially 
assumed parameters. Application of the proposed 
method to a two-dimensional, minimum -time, Earth- 
Mars transfer example reveals a significant reduc- 
tion in computer time required for convergence, when 
compared to the previously proposed quasilineariza- 
tion methods. 
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A MODIFIED QUASILINEARIZATION CONCEPT FOR SOLVING 

THE NONLINEAR TWO-POINT BOUNDARY VALUE PROBLEM 

By Jay M. Lewallen 
Manned Spacecraft Center 


SUMMARY 


A new quasilinearization method, referred to as the Modified Quasilinearization 
Method, is proposed for numerically solving the nonlinear two -point boundary value 
problem. Application of the proposed method to a two-dimensional, minimum-time, 
Earth -Mars transfer reveals a significant reduction in computer time required for con- 
vergence, when compared to previously proposed quasilinearization methods. More- 
over, the iterative method suggested substantially reduces the sensitivity of the 
convergence characteristics of the method to required initially assumed parameters. 

For the case shown, the proposed method required 69 percent less computer time 
than the Generalized Newton-Raphson Method. Also, the suggested iteration method 
increases the convergence envelope, representing the initially assumed values, by 
approximately 350 percent. 


INTRODUCTION 


The optimization of spacecraft trajectories has been of considerable interest for 
a number of years, and significant progress has been made in building a capability for 
solving very complex trajectory problems. One class of trajectory optimization prob- 
lems occurs when it is desired to determine the history of the variables which are ca- 
pable of controlling the spacecraft state in such a manner that certain specified terminal 
constraints are satisfied while some index of performance is extremized. The calculus 
of variations is a classical tool for solving such problems, and with its use the optimi- 
zation problem may be reduced to a two -point boundary value problem. Methods for 
solving the variational two-point boundary value problem are designated indirect meth- 
ods. The convergence characteristics of these methods are extremely sensitive to ini- 
tial Lagrange multiplier values. However, if convergence occurs it is rapid, and in 
some cases the convergence is quadratic. The indirect methods have been successfully 
investigated by Jurovics and McIntyre (ref. 1); Breakwell, Speyer, and Bryson (ref. 2); 
Jazwinski (ref. 3); and McGill and Kenneth (ref. 4). 

The present investigation makes an extension to the quasilinearization concepts 
as presented by McGill and Kenneth (ref. 4). This extension places the quasilineari- 
zation approach in a more competitive position with the other indirect methods; first, 



by allowing the specification of the terminal boundary in terms of a general function of 
the problem variables and time rather than just specific values of the variables them- 
selves, and second, by making the terminal time determination an integral part of the 
iteration process. This time determination is made without requiring any additional 
terms to be added to the existing differential equations, and no additional differential 
equations are needed. 

In addition to the previously stated extensions, an iteration procedure is discussed 
which substantially reduces the convergence sensitivity to initially assumed parameters, 
thus making the method more competitive with the direct methods. This sensitivity re- 
duction is made by requesting only a percentage of the terminal constraint dissatisfac- 
tion to be corrected on a given iteration. 

The Modified Quasilinearization Method has been successfully implemented by 
solving a low constant thrust, two-dimensional, minimum-time, Earth-Mars transfer. 
The method is currently being applied to a rendezvous problem having time- dependent 
terminal constraints. 


SYMBOLS 

A 2n by 2n matrix of partial derivatives defined in equation (2) 
a additional state variable defined in equation (10) 

B 2n vector of partial derivatives defined in equation (2) 

C 2n - p vector of constant corrections 
c iteration factor 

F 2n vector function representing the equations of motion and Euler- Lagrange 
equations 

GM gravitational constant of the sun 
g p vector of general initial boundary conditions 

h 2n + 1 - p vector of general terminal boundary conditions 

m spacecraft mass 

N summation index 

p number of initially specified conditions 
r radial position 

s new independent variable defined in equation (10) 
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T thrust magnitude 

t scalar independent variable time 

u radial velocity 

v tangential velocity 

w 2n vector of dependent variables of the particular differential equation defined 
in equation (7) 

Y 2n x 2n - p matrix of solutions of equation (6) 

y 2n vector of dependent variables of the homogeneous differential equation defined 
in equation (6) 

z 2n vector composed of n state variables and n Euler -Lagrange variables 

j3 control variable, thrust direction relative to the local horizontal 

0 angular position 

A n vector of Euler -Lagrange variables 

p scalar metric defined in equation (9) 

Subscripts: 
f terminal time 

1 index 

n nth trajectory 

0 initial time 

Superscripts: 

1 index 

( ) assumed value 
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Operators: 

( ) total derivative with respect to time 
( )' total derivative with respect to s 
6( ) variation operator 


QUASILINEARIZATION CONCEPTS 


The indirect methods proposed by Jurovics and McIntyre (ref. 1); Breakwell, 
Speyer, and Bryson (ref. 2); and Jazwinski (ref. 3) involve the integration of a set of 
nonlinear differential equations of motion. The coefficients for the linear adjoint or 
perturbation differential equations are formed with the variables generated by these 
nonlinear equations. The quasilinearization approach can be formulated by linearizing 
the differential equations of motion, then using the adjoint or perturbation functions in 
the same general manner as before. The coefficients used to generate a new nominal 
trajectory can be formed from the solution that corresponds to the previous nominal 
trajectory. This quasilinearization concept involves the solution of this sequence of 
linear differential equations, the solution of which converges, under appropriate condi- 
tions, to the solution of the desired nonlinear problem. Since the equations are linear, 
the terminal constraints can be satisfied on every iteration, if desired. However, the 
classical optimality condition is not satisfied until convergence has occurred; and even 
though the end points of the trajectory are satisfied, some care must be taken to insure 
that the trajectory shape is correct. Another characteristic of the quasilinearization 
techniques is that an initially assumed solution is required. However, if a reasonable 
estimate of this solution cannot be made, a starting solution (derived from integrating 
the nonlinear differential equations) may be adequate to result in convergence. This 
requires only that the initial values of the unknown variables be assumed, rather than 
the complete solution. 

The problem is formulated in terms of an ordinary, first-order, nonlinear, vec- 
tor differential equation 


z = F(z, t) 


( 1 ) 


where z is a 2n vector composed of n state variables and n Euler -Lagrange vari- 
ables, and t is the independent variable time. This equation may be expanded about 
the previous nominal trajectory, designated as the nth trajectory, and by ignoring the 
nonlinear terms yields 


z = Az + B 


( 2 ) 
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where 


A = f|H (3) 

L3zJ n 

and 

B = z - Az (4) 

n n 


These terms are evaluated from the previous nominal trajectory. 

Suppose that p of the initial values of z are specified, that is, 


z. 

1 



i = 1, . . ., p 


(5) 


This implies that 2n - p initial values of z must be assumed along with an initial 
time t . The homogeneous part of equation (2) may be expressed as 


y = Ay 


( 6 ) 


and hence is similar to the perturbation equation. Equation (6) may be integrated for- 
ward 2n - p times with each successive starting vector consisting of all zero elements 
except for the element that corresponds to one of the unknown initial conditions. This 
element is set equal to unity. This procedure leads to a 2n x 2n - p matrix solution 

Y ('-‘o)' 

The nonhomogeneous solution to equation (2) may be obtained as a solution to 


w = Aw + B (7) 


which generates a particular solution when integrated forward with the p known initial 
conditions and 2n - p assumed conditions. Now the general solution of the linear 
system, equation (2), becomes 


z(t) = Y(t,t Q )C + w(t) 


( 8 ) 
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where C is a 2n - p vector of constants and w(t) is the 2n vector of nonhomogeneous 
solutions to equation (7). 

Since 2n + 1 - p conditions on the terminal value of z must be specified for a 
variable final-time problem, any 2n - p of these conditions may be selected, and the 
appropriate 2n - p members of equation (8) may be evaluated at the assumed terminal 
time. Then these equations are solved for the 2n - p constants C, which are correc- 
tions used to update the assumed initial conditions for the next iteration. 

This procedure continues until a metric (which represents the maximum distance, 
over the coijiplete independent variable range, between successive nominal trajectories) 
becomes less than some preselected value. This metric is defined by 


N 

p = ^2 max 
i=l t 


z 1 1 (t) - z *(t) 
n+1 ' ' n ' ’ 


(9) 


where n refers to the nth iteration. As this metric decreases, the optimal trajectory 
shape is converged upon for the assumed value of the terminal time. The one remain- 
ing unused terminal condition is used in a conventional scalar application of the Newton- 
Raphson iteration technique to produce a more accurate determination of terminal time. 
The above procedures are repeated until the time corrections are smaller than some 
preselected value. 

The major objections to the Generalized Newton- Raphson Method are that the ter- 
minal time determination is very time-consuming (especially when a large error is 
made in the assumed terminal time) and the initial and terminal conditions must be sim- 
ply specific values of the variables z, rather than general functions of these variables. 
The first objection has been partially removed by Long (ref. 5). 

The method proposed by Long, designated as the Modified Generalized Newton- 
Raphson Method, involves a change of the independent variable t = as where a is a 
constant and is considered a new state variable, and s is a new independent variable 
having values 0 ^ s ^ 1. Now, the differential equation, equation (1), may be written 


z’ 


dz 

ds 


= aF(z, as) 


( 10 ) 


The constant a is considered a new state variable and an additional differential equa- 
tion a’ = 0 may be added, but this is clearly not necessary since its solution is trivial 
The value of a is initially assumed and corrected on each iteration and can be seen to 
be the desired terminal time. 


6 



The determination of terminal time now becomes an integral part of the iterative 
procedure; and its separate consideration, as required by the Generalized Newton- 
Raphson Method, is not necessary. However, a relatively complex term that corre- 
sponds to the new state variable a must be added to each differential equation. 
Moreover, the linear system from which the corrections are obtained is increased 
since the value of a must be iteratively determined. 


MODIFIED QUASILINEARIZATION METHOD 


The method proposed here, the Modified Quasilinearization Method, uses the 
quasilinearization concepts previously outlined, and removes both of the stated objec- 
tions. The manner in which the terminal time is determined proves more satisfactory 
than in any other known quasilinearization method, and the generality of the terminal 
constraints makes the method more competitive with the adjoint and perturbation tech- 
niques proposed by Breakwell, Speyer, and Bryson (ref. 2), and by Jazwinski (ref. 3). 

Equation (8) may be rewritten and evaluated at the terminal time to yield 

Y (V *o) C = 2 (‘f) ' w (*f) (11) 


The right-hand side of this equation is the difference between the desired terminal val- 
ue of z and the linear calculation of the terminal value of w. This difference is inter- 
preted as the variation of zftA and is expressed as 6zftA Now, if both sides of 


equation (11) are premultiplied by 



the resulting expression becomes 



( 12 ) 


where 


! f 

of terminal 

self. Since 


is a 2n + 1 - p x 2n matrix describing the partial change of a general set 

boundary conditions h (z^, t^ j to a change in the terminal values of z^J it- 
the right-hand side of equation (12) is 6h^J, a first-order expansion 


dh = 6h + h dtj 


(13) 
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may be introduced to yield 


dh 


= .t Y (V‘o> 


C + h dt 


(14) 


where dh is a 2n + 1 - p vector of terminal constraint dissatisfaction and dtj is the 
time correction to be made for the next iteration. 

Equation (14) is analogous to the linear algebraic equation derived by Jazwinski 
(ref. 3), the difference being that in the present case the nonlinear differential equations 
of motion are linearized. When the optimization problem is reduced to a nonlinear two- 
point boundary value problem, p becomes equal to n, and the implementation of the 
two methods is similar. A detailed presentation of the numerical procedure described 
earlier is made by Lewallen (ref. 6). 

The iteration philosophy is that only a percentage of the terminal dissatisfaction 
be requested for correction on any given iteration, that is, dh = -ch where 0 ^ c ^ 1. 

It is expected that if a 100-percent correction is requested (normal iteration scheme) 
in cases where the linear representation is poor, the sequence of linear solutions will 
diverge. The less severe request of only a percentage correction is applied initially, 
and successively larger percentages are requested after each convergent iteration. 

Upon each divergent iteration, if any, the percentage correction may be reduced. 


APP LIC ATION 


The system model to be considered is that of minimizing the transfer time of a 
constant low-thrust rocket traveling between the circular and coplanar orbits of Earth 
and Mars. The associated differential and algebraic equations are outlined in the ap- 
pendix. The single control variable is the thrust direction relative to the local horizon- 
tal. 


Figure 1 illustrates how the metric p is reduced as a function of IBM 7094 com- 
putation time for a typical example using the Generalized Newton- Raphson Method, 
Modified Generalized Newton- Raphson Method, and Modified Quasilinearization Method. 
The cases shown are those for initial errors in the two Euler- Lagrange variables A^ 

and X 20 and in the terminal time t f of -10, -10 and 20 percent, respectively. The 

third Euler -Lagrange variable is initially set to negative unity for scaling purposes, 
and A^ is easily determined to be zero. The initial solution for these cases are ob- 
tained from a linear approximation of the solutions. 

The Generalized Newton-Raphson Method, as shown in figure 1(a), requires con- 
vergence to an acceptable metric for an assumed terminal time. A time iteration is 
then made, and the previously reduced metric is degraded to some extent. This 
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Metric p 



0 10 20 30 40 50 60 

Computation time, seconds 



(b) Modified Generalized Newton- 
fa) Generalized Newton- Raphson Raphson Method and Modified 

Method. Quasilinearization Method. 

Figure 1. - Metric p as a function of computation time using the linear initial 
solution, normal iteration scheme, and errors = -10 percent, 

SXgo = -10 percent, and 6t^ = 20 percent. 

degradation is indicated on the plots by dashed lines. The Modified Quasilinearization 
Method, for the case shown, is clearly more satisfactory in terms of required compu- 
tational time, requiring 69 percent less than the time required for the Generalized 
Newton- Raphson Method and 28 percent less than the time required for the Modified 
Generalized Newton-Raphson Method. It should be noted that in the terminal phases of 
convergence, the convergence rate is almost quadratic. 

Many cases have been investigated in an effort to determine how sensitive the 
method is to poor initial assumption for the Euler -Lagrange variables and terminal 
time. These numerical results are best illustrated by building envelopes of conver- 
gence, the boundary of which represents the last convergent trial. The errors that are 
discussed are the percentage errors or deviations from the values that result in an op- 
timal solution. 

Figure 2 illustrates the convergence envelopes for the cases where the optimal 
terminal time is initially assumed, and the initial solution is obtained by integrating the 
nonlinear differential equations of motion. A 100-percent correction in the terminal 
constraints is requested on every iteration in figure 2(a). In figure 2(b), a 50-percent 
correction is requested initially; and this factor is either increased or decreased by 
10 percent on each successive iteration, depending on the convergence or divergence, 
respectively, of the previous trajectory. It is seen that when the initial value of the 
iteration factor is reduced from 100 to 50 percent, the convergence envelope size is 
increased by approximately 350 percent. The primary significance of this result is 
that, by reducing the initial iteration factor, the chances of obtaining a convergent so- 
lution with only one computer run are enhanced greatly. Further increase in the 
convergence envelope size is realized by decreasing the initial value of the iteration 
factor. 
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6A IQ, percent 



-100 ' -50 0 50 

6Xjq, percent 


(a) Initial iteration factor, (b) Initial iteration factor, 

100 percent. 50 percent. 

Note: The numbers indicate the iterations required for convergence. 

Figure 2. - Convergence envelope for the Modified Quasilinearization Method us- 
ing a nonlinear initial solution, terminal time error of zero percent, and initial 
iteration factors of 100 and 50 percent. 


It might be speculated that, since the optimal terminal time is used as an initial 
condition for the first iteration, an easy case has been selected for illustration. Actu- 
ally, cases were examined where terminal time errors were -20 and +20 percent. 

The convergence envelopes for -20 percent terminal time error were smaller than 
the envelopes shown in figure 2, but the envelopes for +20 percent terminal time 
error were larger as indicated in more detail by Lewallen (ref. 6). Reference 6 also 
shows that the Modified Quasilinearization Method compares very favorably with the 
method proposed by Jazwinski (ref. 3). 


CONCLUSIONS 


A theoretical extension of the quasilinearization concept, as applied to the nonlin- 
ear two-point boundary value problem, is possible which allows the terminal boundary 
to be specified as a general function of the problem variables. The numerical results 
indicate that for variable final time problems the suggested method of terminal time 
determination represents a significant reduction of computational time relative to pre- 
viously published methods. Moreover, the proposed iteration philosophy is such that the 
convergence sensitivity of the method to the initially assumed parameters is greatly re- 
duced, thus allowing convergence to occur for many cases which otherwise would have 
diverged. 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, April 7, 1967 
039-00-00-82-72 
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APPENDIX 


The following application is formulated to illustrate the procedure discussed in 
the text. The problem is formulated in terms of the nonlinear differential equations of 
motion for the two-dimensional, Earth-Mars transfer under the influence of a constant 
low thrust 



where u and v are the radial and tangential velocities, respectively, and r and 9 
are the radial and angular displacements, respectively. The angle /3 is the angle the 
thrust vector T makes with the local horizontal, and m is the spacecraft mass. 

The linear Euler- Lagrange differential equations must also be satisfied for a 
trajectory optimization problem 



11 


The classical optimality condition is used to eliminate the control variable terms 
in the differential equations of motion. This relation yields 


sin /3 = 




cos /3 = 


-V^ 


(17) 


The specified initial conditions are 


g l - u ( t o) - u o = 0 
g 2 ■ v (*o) - v o - 0 
g 3 ' r (‘o) - r o = 0 
*4 ' 9 (*o) - e o ' 0 


(18) 


where t is specified. 


The specified terminal conditions are 


h i - u (‘f) - u t - 0 


h 2 ' v l f - v f * 0 


h 3 * r ( l t) - r f * 0 


( 19 ) 


If it is desired to minimize terminal time, two additional terminal conditions are de- 
rived from the transversality conditions. 


h 4 “ X 4f = 0 


h 5 = 1 + X 1 F 1 + X 2 F 2 + X 3 F 3 + X 4 F 4 = 0 


( 20 ) 
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The Lagrange multipliers are scaled by requiring that = -1; hence, the last 
terminal boundary condition hg is ignored. The boundary conditions become 



For the solution of the eight differential equations, 10 boundary conditions must 
be known. If the initial time is assumed, the above nine relations are adequate for 
solution. 

The nonhomogeneous, linear, vector differential equation z = Az + B is com- 
posed of n = 4 linearized differential equations of motion (with the control terms 
eliminated by use of the optimality condition) and n = 4 linearized Euler- Lagrange 
differential equations. These equations are 
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( 29 ) 






letting GM = gravitational constant for the Earth, T = spacecraft thrust, and 
m = m Q - mt = spacecraft mass. 
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These nonhomogeneous linear equations are integrated from t to t, with the 

starting vector z (t Q ) = u v r 6 X^ X^ X 3 xj , where the bar indicates an 

n+l t Q 

assumed value. The homogeneous linear equations (same as above except without 
the , i = 1, . . . , 2n terms) are integrated from t Q to t^ with starting vectors 
n 


y l 



~0~ 


~0~ 


"o' 

0 


0 


0 

0 


0 


0 

0 

/, \ 

0 


0 


y 2 y 


y 3 ( t 0 ) 


1 

Z ' 0) n+l 

0 

6 ''°' n+l 

0 

0 


1 


0 

0 


0 


0 

_0_ 


_0_ 


_1_ 


(31) 


When the assumed terminal time t^ is reached, the algebraic equations that 


must be solved for the corrections to be applied to the initially assumed parameters 


10 ’ 20 ’ 


X. n and t , are 
40 f 


1 

H* 

O 

1 


~ y ll 

y 12 

y 13 


-1 

~du “ 

6X 20 


y 21 

y 22 

y 23 

v f 


dv 

6X 40 


y 31 

y 32 

y 33 

r f 


dr 

1 

O 

1 


_ y 81 

y 82 

00 

CO 

X 4f_ 


dx 4 


(32) 


where the elements in the matrix are evaluated at t^. These corrections are applied, 
and a new nominal trajectory is integrated using z = Az + B. 
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The numerical implementation used the following values on an Earth-Mars trans- 


fer example. 

1 2 

Astronomical unit (AU), meters 0. 14959870 x 10 

Orbital radius of Earth, r E> AU 0. 10000000 x 10 1 

Orbital radius of Mars, r M , AU 0. 15236790 x 10 1 

3/3 21 

Gravitational constant of the Sun, GMg, m / sec .... 0. 13271504 x 10 

3 

Initial spacecraft mass, m Q , kilograms 0. 67978852 x 10 

Thrust, T, newtons 0. 40312370 x 10* 

Mass rate, m, kg/sec 0. 10123858 x 10”^ 
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